library(ggplot2)
  library(reshape2)
  library(dplyr)
  library(tidyr)
  library(GGally)
  library(grid)
    "%&%" = function(a,b) paste(a,b,sep="")
  my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
  out.dir <- '~/GitHub/GenArch/GenArchPaper/OTD_enrichment/'

Pull h2 estimate CI>0 genes for each tissue & rm .\d+ from ensid

for(tistype in c("Tissue-Specific","Tissue-Wide")){
  tsfile <- my.dir %&% 'GTEx_' %&% tistype %&% '_local_h2_se_geneinfo_no_description.txt'
  ts <- read.table(tsfile,sep='\t',header=T)
  tislist <- scan(my.dir %&% 'GTEx_nine_tissues','c')
  tislist <- gsub("-","",tislist)
  tislist <- c("CrossTissue",tislist)
  h2tislist <- "h2." %&% tislist
  setislist <- "se." %&% tislist
  geneinfo <- ts[,1:3]
  for(i in seq_along(tislist)){
    h2tis <- h2tislist[i]
    setis <- setislist[i]
    data <- ts %>% select(AssociatedGeneName, EnsemblGeneID, matches(h2tis), matches(setis)) %>% mutate(ensid=substr(EnsemblGeneID,1,15))
    cidata <- data[data[,3]-data[,4]*2>0,] ##pull genes with non-zero confidence intervals
    print(tislist[i] %&% " " %&% tistype %&% ": " %&% dim(cidata)[1] %&% " non-zero h2 CI genes")
    write.table(cidata, file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_non-zeroCIgenes_info.txt",quote=FALSE,row.names = FALSE)
    write.table(cidata[,5], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_non-zeroCIgenes_ensid_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
    write.table(cidata[,1], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_non-zeroCIgenes_gene_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
    
    for(thresh in c(0,0.01,0.05,0.1)){
      h2data <- data[data[,3]>thresh,] ##pull genes with h2 > thresh
    print(tislist[i] %&% " " %&% tistype %&% ": " %&% dim(h2data)[1] %&% " h2 > " %&% thresh %&% " genes")
    write.table(h2data, file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_h2_" %&% thresh %&% "genes_info.txt",quote=FALSE,row.names = FALSE)
    write.table(h2data[,5], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_h2_" %&% thresh %&% "genes_ensid_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
    write.table(h2data[,1], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_h2_" %&% thresh %&% "genes_gene_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
    } 
  }
}
## [1] "CrossTissue Tissue-Specific: 2938 non-zero h2 CI genes"
## [1] "CrossTissue Tissue-Specific: 17007 h2 > 0 genes"
## [1] "CrossTissue Tissue-Specific: 10282 h2 > 0.01 genes"
## [1] "CrossTissue Tissue-Specific: 5042 h2 > 0.05 genes"
## [1] "CrossTissue Tissue-Specific: 2715 h2 > 0.1 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 207 non-zero h2 CI genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 17007 h2 > 0 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 6908 h2 > 0.01 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 2463 h2 > 0.05 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 650 h2 > 0.1 genes"
## [1] "ArteryTibial Tissue-Specific: 306 non-zero h2 CI genes"
## [1] "ArteryTibial Tissue-Specific: 17007 h2 > 0 genes"
## [1] "ArteryTibial Tissue-Specific: 7779 h2 > 0.01 genes"
## [1] "ArteryTibial Tissue-Specific: 3126 h2 > 0.05 genes"
## [1] "ArteryTibial Tissue-Specific: 892 h2 > 0.1 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 298 non-zero h2 CI genes"
## [1] "HeartLeftVentricle Tissue-Specific: 17007 h2 > 0 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 7796 h2 > 0.01 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 4012 h2 > 0.05 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 1614 h2 > 0.1 genes"
## [1] "Lung Tissue-Specific: 223 non-zero h2 CI genes"
## [1] "Lung Tissue-Specific: 17007 h2 > 0 genes"
## [1] "Lung Tissue-Specific: 7143 h2 > 0.01 genes"
## [1] "Lung Tissue-Specific: 2738 h2 > 0.05 genes"
## [1] "Lung Tissue-Specific: 797 h2 > 0.1 genes"
## [1] "MuscleSkeletal Tissue-Specific: 434 non-zero h2 CI genes"
## [1] "MuscleSkeletal Tissue-Specific: 17007 h2 > 0 genes"
## [1] "MuscleSkeletal Tissue-Specific: 7371 h2 > 0.01 genes"
## [1] "MuscleSkeletal Tissue-Specific: 2356 h2 > 0.05 genes"
## [1] "MuscleSkeletal Tissue-Specific: 622 h2 > 0.1 genes"
## [1] "NerveTibial Tissue-Specific: 379 non-zero h2 CI genes"
## [1] "NerveTibial Tissue-Specific: 17006 h2 > 0 genes"
## [1] "NerveTibial Tissue-Specific: 7797 h2 > 0.01 genes"
## [1] "NerveTibial Tissue-Specific: 3491 h2 > 0.05 genes"
## [1] "NerveTibial Tissue-Specific: 1216 h2 > 0.1 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 449 non-zero h2 CI genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 17007 h2 > 0 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 8604 h2 > 0.01 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 3589 h2 > 0.05 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 1097 h2 > 0.1 genes"
## [1] "Thyroid Tissue-Specific: 442 non-zero h2 CI genes"
## [1] "Thyroid Tissue-Specific: 17007 h2 > 0 genes"
## [1] "Thyroid Tissue-Specific: 7757 h2 > 0.01 genes"
## [1] "Thyroid Tissue-Specific: 3141 h2 > 0.05 genes"
## [1] "Thyroid Tissue-Specific: 1010 h2 > 0.1 genes"
## [1] "WholeBlood Tissue-Specific: 490 non-zero h2 CI genes"
## [1] "WholeBlood Tissue-Specific: 17007 h2 > 0 genes"
## [1] "WholeBlood Tissue-Specific: 8504 h2 > 0.01 genes"
## [1] "WholeBlood Tissue-Specific: 2918 h2 > 0.05 genes"
## [1] "WholeBlood Tissue-Specific: 819 h2 > 0.1 genes"
## [1] "CrossTissue Tissue-Wide: 2938 non-zero h2 CI genes"
## [1] "CrossTissue Tissue-Wide: 17007 h2 > 0 genes"
## [1] "CrossTissue Tissue-Wide: 10282 h2 > 0.01 genes"
## [1] "CrossTissue Tissue-Wide: 5042 h2 > 0.05 genes"
## [1] "CrossTissue Tissue-Wide: 2715 h2 > 0.1 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 1127 non-zero h2 CI genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 17007 h2 > 0 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 6236 h2 > 0.01 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 3407 h2 > 0.05 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 1736 h2 > 0.1 genes"
## [1] "ArteryTibial Tissue-Wide: 1171 non-zero h2 CI genes"
## [1] "ArteryTibial Tissue-Wide: 17007 h2 > 0 genes"
## [1] "ArteryTibial Tissue-Wide: 6314 h2 > 0.01 genes"
## [1] "ArteryTibial Tissue-Wide: 3546 h2 > 0.05 genes"
## [1] "ArteryTibial Tissue-Wide: 1889 h2 > 0.1 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 749 non-zero h2 CI genes"
## [1] "HeartLeftVentricle Tissue-Wide: 17007 h2 > 0 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 6352 h2 > 0.01 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 3865 h2 > 0.05 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 2179 h2 > 0.1 genes"
## [1] "Lung Tissue-Wide: 929 non-zero h2 CI genes"
## [1] "Lung Tissue-Wide: 17007 h2 > 0 genes"
## [1] "Lung Tissue-Wide: 6331 h2 > 0.01 genes"
## [1] "Lung Tissue-Wide: 3387 h2 > 0.05 genes"
## [1] "Lung Tissue-Wide: 1693 h2 > 0.1 genes"
## [1] "MuscleSkeletal Tissue-Wide: 1063 non-zero h2 CI genes"
## [1] "MuscleSkeletal Tissue-Wide: 17007 h2 > 0 genes"
## [1] "MuscleSkeletal Tissue-Wide: 6049 h2 > 0.01 genes"
## [1] "MuscleSkeletal Tissue-Wide: 2814 h2 > 0.05 genes"
## [1] "MuscleSkeletal Tissue-Wide: 1331 h2 > 0.1 genes"
## [1] "NerveTibial Tissue-Wide: 1466 non-zero h2 CI genes"
## [1] "NerveTibial Tissue-Wide: 17006 h2 > 0 genes"
## [1] "NerveTibial Tissue-Wide: 6594 h2 > 0.01 genes"
## [1] "NerveTibial Tissue-Wide: 4057 h2 > 0.05 genes"
## [1] "NerveTibial Tissue-Wide: 2399 h2 > 0.1 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 1198 non-zero h2 CI genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 17007 h2 > 0 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 6504 h2 > 0.01 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 3530 h2 > 0.05 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 1809 h2 > 0.1 genes"
## [1] "Thyroid Tissue-Wide: 1327 non-zero h2 CI genes"
## [1] "Thyroid Tissue-Wide: 17006 h2 > 0 genes"
## [1] "Thyroid Tissue-Wide: 6590 h2 > 0.01 genes"
## [1] "Thyroid Tissue-Wide: 3758 h2 > 0.05 genes"
## [1] "Thyroid Tissue-Wide: 2082 h2 > 0.1 genes"
## [1] "WholeBlood Tissue-Wide: 945 non-zero h2 CI genes"
## [1] "WholeBlood Tissue-Wide: 17007 h2 > 0 genes"
## [1] "WholeBlood Tissue-Wide: 5426 h2 > 0.01 genes"
## [1] "WholeBlood Tissue-Wide: 2586 h2 > 0.05 genes"
## [1] "WholeBlood Tissue-Wide: 1261 h2 > 0.1 genes"
write(data$ensid, file=out.dir %&% "full_tested_ensid_list.txt",ncolumns=1)

Define known genes for 7 WTCCC diseases based on the GWAS catalog and make list of ALL genes in catalog

  #catalog used:
  gwasfile <- my.dir %&% 'gwas_catalog_v1.0-downloaded_2015-06-02.tsv'
#  gwas <- data.frame(read.table(gwasfile,header=TRUE,sep='\t',quote="",comment.char="",as.is=TRUE)) #for reference
  runPERL <- "perl " %&% my.dir %&% "24_define_disease_genes.pl " %&% gwasfile
  system(runPERL)

Disease gene enrichment in top h2 genes by tissue

lci<-function(x) quantile(x, c(.025, 0.975),na.rm=T)[[1]]
uci<-function(x) quantile(x, c(.025, 0.975),na.rm=T)[[2]]
##is testvec signficanty enriched for setvec? make n samples of size(testvec) from fullvec and count overlap
enrichment <- function(setvec, testvec, fullvec, nperms = 1000){
  obs <- length(testvec[testvec %in% setvec])
  counts <- vector()
  for(i in 1:nperms){
    rantest <- base::sample(fullvec,length(testvec),replace=FALSE)
    cnt <- length(rantest[rantest %in% setvec])
    counts <- c(counts,cnt)
  }
  empp <- mean(counts>obs)
  meanc <- mean(counts)
  lc <- lci(counts)
  uc <- uci(counts)
  return(c(obs,meanc,lc,uc,empp))
}


  
set.seed(42)
fullgenelist <- as.character(ts$AssociatedGeneName)

dislist <- c("BD","CAD","HT","T1D","T2D","CD","RA","ALL")
nperms <- 10000 

for(thresh in c("non-zeroCI","h2_0.1","h2_0.05","h2_0.01")){
  results <- data.frame(type=character(0),dis=character(0),tis=character(0),obsOverlap=double(0),meanOverlap=double(0),lCI=double(0),uCI=double(0),empPval=double(0))
  for(tistype in c("Tissue-Specific","Tissue-Wide")){
    for(i in seq_along(dislist)){
      dis <- dislist[i]
      disgenes <- scan(out.dir %&% "gwas." %&% dis %&% ".tsv","character")
      for(j in seq_along(tislist)){
        tis <- tislist[j]
        tisgenes <- scan(out.dir %&% tis %&% '_' %&% tistype %&% '_' %&% thresh %&% 'genes_gene_list.txt','character')
        res <- enrichment(disgenes,tisgenes,fullgenelist,nperms=nperms)
        resvec <- data.frame(type=tistype,dis=dis,tis=tis,obsOverlap=res[1],meanOverlap=res[2],lCI=res[3],uCI=res[4],empPval=res[5])
        results <- rbind(results,resvec)
      }
    }
  }
  sortres <- arrange(results,empPval)
  write.table(sortres,file=out.dir %&% "GWAS_catalog_disease_gene5e-8_enrichment_in_" %&% thresh %&% "_genes_by_tissue_" %&% nperms %&% "perms.txt",row.names=FALSE,quote=FALSE)
  print("GWAS catalog disease gene enrichment in " %&% thresh %&% " genes by tissue " %&% nperms %&% " perms:")
  print(head(sortres,n=20L))
}
## [1] "GWAS catalog disease gene enrichment in non-zeroCI genes by tissue 10000 perms:"
##               type dis                      tis obsOverlap meanOverlap lCI
## 1  Tissue-Specific ALL                     Lung         82     56.2072  44
## 2      Tissue-Wide ALL               WholeBlood        290    238.4374 213
## 3  Tissue-Specific ALL                  Thyroid        140    111.5018  94
## 4  Tissue-Specific  CD           MuscleSkeletal         17      8.2607   3
## 5  Tissue-Specific T2D                     Lung          5      1.1964   0
## 6  Tissue-Specific  HT                     Lung          2      0.2605   0
## 7  Tissue-Specific T1D               WholeBlood          6      2.1424   0
## 8  Tissue-Specific CAD       HeartLeftVentricle          4      1.3001   0
## 9  Tissue-Specific CAD           MuscleSkeletal          5      1.9170   0
## 10 Tissue-Specific ALL           MuscleSkeletal        128    109.6354  92
## 11 Tissue-Specific ALL       HeartLeftVentricle         90     75.1409  61
## 12     Tissue-Wide ALL                     Lung        260    234.2492 209
## 13 Tissue-Specific T2D SkinSunExposed(Lowerleg)          5      2.4301   0
## 14 Tissue-Specific T1D       HeartLeftVentricle          3      1.2897   0
## 15 Tissue-Specific  RA       HeartLeftVentricle          5      2.4808   0
## 16 Tissue-Specific ALL               WholeBlood        139    123.5206 105
## 17 Tissue-Specific  HT       HeartLeftVentricle          1      0.3530   0
## 18 Tissue-Specific  RA               WholeBlood          7      4.0554   1
## 19 Tissue-Specific T2D               WholeBlood          5      2.6384   0
## 20 Tissue-Specific ALL             ArteryTibial         88     77.2368  63
##        uCI empPval
## 1   69.000  0.0000
## 2  264.000  0.0003
## 3  129.025  0.0009
## 4   14.000  0.0013
## 5    4.000  0.0015
## 6    2.000  0.0020
## 7    5.000  0.0042
## 8    4.000  0.0098
## 9    5.000  0.0111
## 10 128.000  0.0206
## 11  90.000  0.0215
## 12 260.000  0.0223
## 13   6.000  0.0351
## 14   4.000  0.0390
## 15   6.000  0.0413
## 16 142.000  0.0444
## 17   2.000  0.0484
## 18   8.000  0.0499
## 19   6.000  0.0508
## 20  92.000  0.0687
## [1] "GWAS catalog disease gene enrichment in h2_0.1 genes by tissue 10000 perms:"
##               type dis                      tis obsOverlap meanOverlap
## 1      Tissue-Wide ALL               WholeBlood        366    318.1679
## 2  Tissue-Specific T1D               WholeBlood          8      3.5559
## 3  Tissue-Specific  CD           MuscleSkeletal         20     11.8198
## 4  Tissue-Specific  HT                     Lung          3      0.9179
## 5  Tissue-Specific CAD           MuscleSkeletal          6      2.7539
## 6      Tissue-Wide ALL             ArteryTibial        512    476.7326
## 7  Tissue-Specific ALL                  Thyroid        278    254.8368
## 8  Tissue-Specific  BD              NerveTibial          3      1.3102
## 9  Tissue-Specific  CD       HeartLeftVentricle         39     30.5712
## 10 Tissue-Specific  HT              NerveTibial          3      1.4320
## 11 Tissue-Specific ALL           MuscleSkeletal        174    156.8398
## 12     Tissue-Wide  HT              NerveTibial          5      2.8251
## 13     Tissue-Wide  CD               WholeBlood         31     23.9342
## 14 Tissue-Specific ALL               WholeBlood        225    206.5587
## 15     Tissue-Wide  BD                  Thyroid          4      2.2264
## 16 Tissue-Specific ALL       HeartLeftVentricle        430    407.0677
## 17 Tissue-Specific ALL                     Lung        217    201.2189
## 18     Tissue-Wide ALL                     Lung        449    426.8637
## 19 Tissue-Specific T2D                  Thyroid          8      5.4615
## 20 Tissue-Specific  BD SkinSunExposed(Lowerleg)          2      1.1466
##        lCI uCI empPval
## 1  289.000 347  0.0005
## 2    1.000   7  0.0080
## 3    6.000  19  0.0086
## 4    0.000   3  0.0112
## 5    0.000   6  0.0182
## 6  442.000 512  0.0245
## 7  229.000 282  0.0405
## 8    0.000   4  0.0413
## 9   21.000  41  0.0474
## 10   0.000   4  0.0482
## 11 136.000 178  0.0505
## 12   0.000   6  0.0520
## 13  15.000  33  0.0536
## 14 183.000 230  0.0585
## 15   0.000   5  0.0615
## 16 374.000 440  0.0799
## 17 178.000 225  0.0857
## 18 394.000 460  0.0907
## 19   1.975  10  0.0916
## 20   0.000   3  0.0996
## [1] "GWAS catalog disease gene enrichment in h2_0.05 genes by tissue 10000 perms:"
##               type dis                      tis obsOverlap meanOverlap lCI
## 1  Tissue-Specific  HT                     Lung          7      3.2244   0
## 2      Tissue-Wide  BD      AdiposeSubcutaneous          7      3.6278   1
## 3      Tissue-Wide ALL           MuscleSkeletal        754    709.8752 669
## 4  Tissue-Specific  BD              NerveTibial          7      3.7060   1
## 5      Tissue-Wide ALL               WholeBlood        692    652.2293 613
## 6      Tissue-Wide ALL              NerveTibial       1070   1023.2204 975
## 7      Tissue-Wide ALL             ArteryTibial        938    894.7990 850
## 8  Tissue-Specific T2D                     Lung         20     14.8128   8
## 9      Tissue-Wide ALL                  Thyroid        984    947.6142 904
## 10 Tissue-Specific ALL                  Thyroid        826    792.2935 749
## 11 Tissue-Specific  HT      AdiposeSubcutaneous          5      2.8948   0
## 12 Tissue-Specific  BD SkinSunExposed(Lowerleg)          6      3.7760   1
## 13 Tissue-Specific CAD      AdiposeSubcutaneous         15     10.8874   5
## 14 Tissue-Specific  BD               WholeBlood          5      3.0924   0
## 15 Tissue-Specific  CD             ArteryTibial         69     59.3060  46
## 16     Tissue-Wide  BD                  Thyroid          6      3.9805   1
## 17 Tissue-Specific CAD               WholeBlood         17     12.8854   7
## 18 Tissue-Specific ALL           MuscleSkeletal        621    594.5320 556
## 19     Tissue-Wide T2D             ArteryTibial         24     19.1016  12
## 20 Tissue-Specific  BD                  Thyroid          5      3.3069   0
##         uCI empPval
## 1     7.000  0.0089
## 2     7.000  0.0178
## 3   751.000  0.0183
## 4     7.000  0.0200
## 5   692.000  0.0240
## 6  1071.025  0.0274
## 7   940.000  0.0284
## 8    22.000  0.0559
## 9   993.000  0.0567
## 10  836.000  0.0581
## 11    6.000  0.0613
## 12    7.000  0.0632
## 13   17.000  0.0711
## 14    6.000  0.0716
## 15   73.000  0.0738
## 16    8.000  0.0800
## 17   20.000  0.0820
## 18  633.000  0.0852
## 19   27.000  0.0865
## 20    7.000  0.0976
## [1] "GWAS catalog disease gene enrichment in h2_0.01 genes by tissue 10000 perms:"
##               type dis                      tis obsOverlap meanOverlap
## 1      Tissue-Wide ALL               WholeBlood       1443   1368.5416
## 2      Tissue-Wide  BD              NerveTibial         12      6.9685
## 3      Tissue-Wide ALL             ArteryTibial       1658   1592.5260
## 4      Tissue-Wide  BD                  Thyroid         11      6.9748
## 5  Tissue-Specific ALL              NerveTibial       2028   1966.7564
## 6      Tissue-Wide ALL           MuscleSkeletal       1582   1525.7863
## 7  Tissue-Specific  BD SkinSunExposed(Lowerleg)         13      9.1227
## 8  Tissue-Specific T2D                     Lung         48     38.7017
## 9  Tissue-Specific  BD             ArteryTibial         12      8.2039
## 10     Tissue-Wide ALL              NerveTibial       1715   1662.6956
## 11 Tissue-Specific  HT              NerveTibial         13      9.1965
## 12     Tissue-Wide ALL                  Thyroid       1715   1662.3288
## 13     Tissue-Wide  BD               WholeBlood          9      5.7428
## 14     Tissue-Wide CAD           MuscleSkeletal         34     26.6712
## 15     Tissue-Wide  BD      AdiposeSubcutaneous         10      6.6147
## 16 Tissue-Specific  CD      AdiposeSubcutaneous        145    131.2328
## 17 Tissue-Specific ALL                  Thyroid       2000   1956.6711
## 18     Tissue-Wide T2D      AdiposeSubcutaneous         40     33.7749
## 19 Tissue-Specific ALL      AdiposeSubcutaneous       1780   1741.8881
## 20     Tissue-Wide  BD             ArteryTibial          9      6.7149
##     lCI  uCI empPval
## 1  1317 1421  0.0032
## 2     3   11  0.0033
## 3  1540 1646  0.0077
## 4     3   11  0.0139
## 5  1912 2022  0.0142
## 6  1473 1579  0.0178
## 7     5   13  0.0179
## 8    30   48  0.0189
## 9     4   12  0.0194
## 10 1609 1716  0.0252
## 11    5   14  0.0256
## 12 1608 1716  0.0265
## 13    2   10  0.0303
## 14   19   35  0.0306
## 15    3   11  0.0322
## 16  114  149  0.0524
## 17 1902 2013  0.0620
## 18   25   43  0.0680
## 19 1687 1796  0.0818
## 20    3   11  0.0869

distribution of h2 for disease vs non disease

dislist <- c("BD","CAD","HT","T1D","T2D","CD","RA","ALL")
tislist <- c("CrossTissue","AdiposeSubcutaneous","ArteryTibial","HeartLeftVentricle","Lung","MuscleSkeletal","NerveTibial","SkinSunExposed(Lowerleg)","Thyroid","WholeBlood")
typelist<-c("Tissue-Specific","Tissue-Wide")

for(thresh in c(0.1,0.05,0)){
  for(tistype in typelist){
    for(tis in tislist){
      info <- read.table(out.dir %&% tis %&% '_' %&% tistype %&% '_h2_' %&% thresh %&% 'genes_info.txt',header=T)
      finaldf <- data.frame(AssociatedGeneName=character(0),EnsemblGeneID=character(0),h2=double(0),se=double(0),ensid=character(0),diseaseGene=logical(0L),disease=character(0))
      for(dis in dislist){
        setvec <- scan(out.dir %&% "gwas." %&% dis %&% ".tsv","character")
        disinfo <- info %>% mutate(diseaseGene=(info[,1] %in% setvec),disease=dis)
        colnames(disinfo) <- c("AssociatedGeneName","EnsemblGeneID","h2","se","ensid","diseaseGene","disease")
        finaldf <- rbind(finaldf,disinfo)
      }
      p<-ggplot(finaldf,aes(x=finaldf[,3],fill=diseaseGene,color=diseaseGene)) + facet_wrap(~disease,ncol=2) + geom_density(alpha=0.3) + xlab("h2") + ggtitle(tistype %&% ' ' %&% tis %&% ' h2 > ' %&% thresh)
      print(p)
      p<-ggplot(finaldf,aes(y=finaldf[,3],x=diseaseGene)) + facet_wrap(~disease,ncol=4) + geom_jitter(aes(colour=diseaseGene),alpha=0.3,position = position_jitter(width = .15)) + geom_boxplot(alpha=0,outlier.size=NA) + xlab("diseaseGene") + ylab("h2") + ggtitle(tistype %&% ' ' %&% tis %&% ' h2 > ' %&% thresh) +theme_bw() + theme(legend.position="none") 
      print(p)
    }
  }
}